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ABSTRACT 

We investigate the variability properties of main sequence stars in the first month of Kepler data, using a new astrophysically robust 
systematics correction. We find that the fraction of stars with variability greater than that of the Sun is 60%, which is marginally 
consistent with previous studies, and confirm the trend of increasing variability with decreasing eff'ective temperatures. We define low 
and high variability samples, with a cut corresponding to twice the variability level of the active Sun, and compare the properties 
of the stars belonging to each sample. We show tentative evidence that the more active stars have lower proper motions and may be 
located closer to the galactic plane. We also investigate the frequency content of the variability, finding clear evidence for periodic 
or quasi-periodic behaviour in 16% of stars, and showing that there exist significant diff'erences in the nature of variability between 
spectral types. Of the periodic objects, most A and F stars have short periods (< 2 days) and highly sinusoidal variability, suggestive of 
pulsations, whilst G, K and M stars tend to have longer periods (> 5 days, with a trend towards longer periods at later spectral types) 
and show a mixture of periodic and stochastic variability, indicative of activity. Finally, we use auto-regressive models to characterise 
the stochastic component of the variability, and show that its typical amplitude and time- scale both increase towards later spectral 
types, which we interpret as a corresponding increase in the characteristic size and life-time of active regions. 
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1. Introduction 

The past decade has seen a rapid increase in the rate of dis- 
covery and classification of variable stars, mainly as a result of 
time-domain photometric surveys whose primary goals were to 
search for other phenomena, such as microlensing or planetary 
transits (e.g. the OGLE survey, Udalski et al. 2008 ) . Sources of 
stellar variability are wide ranging, from the potentially large- 
amplitude signatures of eclipses, star spots and pulsations, down 
to the sub-millimagnitude changes induced by granulation. The 
typical precision and time sampling of ground-based surveys 
confines the associated variability studies to 'classical' variable 
stars, with amplitudes of a percent or more. In these surveys, the 
Sun, whose total output never varies by more than around 0.5% 
peak-to-peak ( Frohlich 2011 ), even at the maximum of its activ- 
ity cycle, would appear as a 'quiet' or 'const ant' star . However, 
space-b ased transit surveys such as CoRoT ( Baglin |2003) and 
Kepler ( Borucki et al.|201Q ) are sensitive to 'micro- variability' 
down to and well below the solar level, on timescales ranging 
from minutes to months and, over the entire lifetime of Kepler, 
years. 

Measuring the basic characteristics of the variability (ampli- 
tude, periodicity, etc. . . ) across large samples of stars, and com- 
paring them to stellar parameters such as age, mass and compo- 
sition, is a first step towards a better understanding of the under- 
lying phenomena. Many of the latter are ill-understood, because 
they are related to rotation, convection and magnetism which 
are challenging to model. Variability statistics also have a cru- 
cial impact on exoplanet studies, particularly for radial- velocity 



searches or radial- velocity confirmation of transiting planet can- 
didates (see e.g. |Pont et al.|2011 ). 

Data from the Kepler mission is particularly amenable to sta- 
tistical variability studies because of the instrument's unprece- 
dented photometric precisioi£]and vast field of view (115 deg^). 
The Quarter 1 (Ql) data, which was made public in June 2010, 
has already been studied by |Basri et al.| (2010, hereafter BIO), 
who show that somewhat less than half of the dwarf stars sur- 
veyed by Kepler are more variable than the Sun on timescales 
of up to a month, with the fraction increasing from earlier to 
later spectral types. |Basri et al. ( 2011 hereafter Bl 1) went on to 
demonstrate that periodic variable stars have significantly larger 
amplitudes, as a sample, than aperiodic variables. Finally, Ciardi] 
|et al.| ( |2011[ hereafter Cll) performed a complementary study 
of the same sample using dispersion rather than amplitude as a 
variability statistic, and studying likely dwarfs and giants sepa- 
rately using the st ellar param eters p rovided in the K epler Input 
Catalog (KIC, Brown etal.|201 1 ; Bat alhaet al.| 20TQl). Variability 
statistics have also been determined using the 10 days of com- 
missioning data (QO), with the aim of developing methods to 



characterise and select specific types of variable (Walkowicz & 
|Basri|2010|). 

In Cll, the data were corrected for systematics using the 
Kepler team's Pre -Search Data Conditioning (PDC) method 
( Jenkins et al.||201 0). This inspired us to investigate further the 
apparent bimodality in the variability of dwarf stars observed by 
Kepler, with particular attention to the eflTect of diflferent system- 
atics correction methods. 



* e-mail: Amy . McQuillan@astro . ox . ac . uk 



^ The Kepler mission requires photometric precision of 20 ppm on 
6.5 hr timescales for a Kepler magnitude 12 star, see [Jenkins et al.| 
"20101) for a detailed discussion. 
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The goal of the present paper is to revisit the work of Cll 
and BIO, 11 following the application of a new astrophysically 
robust de-trending method, designed to preserve intrinsic vari- 
ability signals and remove as fully as possible the systematics. 
We quantify trends, previously identified and new, between vari- 
ability characteristics and stellar properties, and investigate the 
nature of the variability in more detail, in order to gain further 
insight into the underlying mechanisms. 

We describe the data and our systematics correction process 
in Section |2] In Section [3] we discuss our choice of variability 
statistics and examine the fractions and physical properties of 
the low and high variability samples, focussing on their periodic 
and stochastic nature in Section |4] We present our conclusions 
and discuss the implications of our findings in Section [5] 



2. Systematics Removal in the Q1 Data 

2.1. The Data 

The Kepler Quarter 1 (Ql) observations took place over -33.5 
days between May 13* and June 15* 2009. 156,097 targets were 
observed during this time, most with a cadence of 29.42 min- 
utes (a small subset were observed with an increased cadence 
of ~1 minute, but only the long cadence data was used in this 
study). The plate scale is 3.98 arcsec per pixel (van Cl eve &| 
|Caldwell|[200 9) and the astrometric precis ion for a single 30 
minute measure is better than 4 milliarcsec ( [Monet et al.|2010| ). 
Known eclip sing binary sy stems were removed based on the list 
compiled by Slawson et a l. (2011), available onlin^ Kepler 
planet candidates ( Borucki et al.| 2011) were also removec^ 
Light curves with discontinuities, detected as a jump greater than 
Act and visually confirmed, were also removed, leaving the final 
sample size at 123030 stars. 

Both the 'raw' and PDC corrected data are publicly avail- 
able through the Kepler mission archiv^ at STScI and also from 
the NASA Star and Exoplanet Database^ (NStED). This study 
uses the STScI August 2011 data release. The stellar proper- 
ties used for classification in this study, namely eflTective tem- 
perature (accurate to 200 K) and surface gravity (accurate to 0.5 
dex), come from the KIC. These parameters were estimated us- 
ing Bayesian posterior probability maximisation to match ob- 
served colors, estimated from Sloan g, r, i, z filters, 2MASS 
JHK, and D51 (510nm), to |Castelli & Kurucz] ( |20Q4| ) stellar at- 
mosphere models. The proper motion values from the KIC are 
taken from a selection of catalog^ where available. Total proper 
motion is listed on NStED as having accuracy of 20 milliarcsec- 
onds per year. 39,000 dwarf stars are listed with non-zero total 
proper motion. The precisions listed in the KIC are the typical 
value for each parameter, although in reality these may vary by a 
small amount between stars of diflTerent magnitude and spectral 
type. For a more detailed discussion of the KIC parameters, see 
[Brown et al.| ( |2QTT] ); |BataTha et al.| ( |20T0| ); |Vemer et aLldMTT] ). 



^ http : //keplerebs . villanova. edu/ 

^ http : //archive. stsci. edu/kepler/planet_candidates .html 

^ http://archive.stsci.edu/kepler 

^ http://nsted.ipac.caltech.edu 
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2.2. Astrophysically Robust Correction (ARC) 

CIO used the PDC corrected flux, whereas B 10,1 1 used the 'raw' 
data, fitting and subtracting a low-order polynomial to remove 
long-term trends before computing their statistics. 

The PDC correction, which is progressively revised and 
updated, uses ancillary engineering data such as temperature, 
pointing and focus variations to remove systematic effects from 
the light curves. The PDC attempts to fit and remove intrinsic 
variability and then apply the systematics correction to the resid- 
uals before re-inserting the stellar signal. This creates a decision 
point where the program must determine whether the variabil- 
ity is real or not, which can lead to real low amplitude and long 
term stellar signal being falsely removed. It can also add high- 
frequency noise ( van Cleve 2010). The polynomial correction of 
BIO, 11 can equally aflTect real long-term variability so we there- 
fore opted to develop a new astrophysically robust correction 
(ARC) for systematics. The ARC will be presented in detail in a 
forthcoming paper (Roberts et al. in prep.), but we give a brief 
outline of it here for completeness. 

Our core assumptions are that i) the trends are systematic 
i.e. they are instrumental in origin and are therefore present, at 
some level, in the majority of light curves, ii) the underlying 
number of significant trends is unknown and the trend profile is 
not pre- specified, iii) the amount of trend present changes from 
light curve to light curve and is unknown and finally iv) although 
many trends may be curves with a long timescale, there should 
be no bias towards smoothness. 

We consider the observed set of light curves, J/, to be com- 
posed of a linear combination of underlying "true" curves, Sf, an 
unknown combination of an unknown number, /, of systematic 
trends, uj, and an observation noise process, 6, 



Si + 



7=1 



Mr + 6; 



(1) 



in which the unknown factors aij represent the amount of trend 
j believed to be present in light curve /. 

To evaluate the number of trends we start by modelling a 
light curve, J„ say, as a linear combination of all other curves 
with factor weights jSi, 



(2) 



Kepler Stellar Classification Program, 
UCAC2, 2MASS and USNO-Bl.O. 



Hipparcos, Tycho-2, 



The inference is achieved using a fully Bayesian linear model, 
based on the computationally efficient variational Bayes method- 
ology and including shrinkage on the /3i parameters. 

If we consider the inferred empirical distribution of all {\l3i\} 
we expect systematic trends to be associated with the highest 
entropy distributions - as true systematic trends are more likely 
to be present in many light curves rather than in just a small 
subset. Significanj^ high-entropy trend components are further 
processed to remove any residual noise, using empirical mode 
decomposition (EMD) ( Huang et al.||1998 ), allowing us to re- 
cover the dominant dynamic of the trend without imposing any 
pre-defined smoothness constraints. The latter enables us to ex- 
tract very low-frequency trends as well as, for example, system- 
atic artefacts due to high frequency vibrations of the satellite. 
After extraction of the dominant dynamics via EMD, the resul- 
tant trend is then removed from all the data by inferring the fac- 
tors Gij in Equation[T] once again using fully Bayesian inference, 

^ We define significance by considering a principal component anal- 
ysis of the highest entropy trends and require the first principal compo- 
nent to explain more that 95% of the spectral radius. 
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Fig. 1. Subplot a) shows example basis trends, inferred from all Ql data; note that the y-axis is in arbitrary units, as these basis 
functions are scaled to support light curves. Subplots b) and c) show examples of de-trending and allow comparison between 
methods: the top trace of each example shows the raw light curve (blue) and the removed trend (red) and the lower subplots show 
the resultant de-trended light curves using the different correction methods. Subplot d) provides detail of a small section of data, 
highlighting the effective satellite vibration artefact removal obtained. 



giving rise to a modified data set which is deflated with respect 
to the trend. The entire process of trend discovery and removal is 
then repeated using the deflated data set until no more significant 
trends are discovered. 

Figure [TJa) shows a typical set of inferred basis trends, ob- 
tained from all light curves over all CCDs; note the two typical 
smooth trends along with the third high frequency vibration ef- 
fect. Plots (b) and (c) show the raw light curves (blue) along with 
the inferred trend component (red) and below the de-trended 
light curves for the different correction methods compared in 
this paper. Plot (d) highlights a typical section of data in which 
the high frequency systematic is prevalent and well-removed by 
the method. These plots demonstrate that the ARC provides the 
most consistently suitable corrections, removing the high fre- 
quency trend in (d), maintaining the true stellar variability in (b) 
and dealing with low frequency smooth systematics effectively 
in (c). Comparatively it can be seen that neither the PDC or B 10 
corrections can perform well in all these areas. 

We compared results obtained when correcting all light 
curves together, with that produced by treating each mod.ouj^ 
separately. The mod.out separated batches are most suitable 
since they use the same CCD and channel, which leads to more 
closely correlated systematics, for example, the reaction motor 
effects are more prevalent on some CCDs than others. 



Each module consists of 2 CCDs, which in turn have 2 outputs each. 



The ARC performs well in the vast majority of cases, remov- 
ing systematic trends without altering intrinsic stellar variability 
signals. One remaining effect that is not currently removed by 
the ARC is that of the variable guide stars, primarily the eclips- 
ing binary. These introduce small variations in some light curves 
but these are not frequently occurring or similar enough to be de- 
tected by the ARC. A method to detect and remove this effect is 
currently being devised and will be included in future work, but 
for this study it may be omitted without significantly reducing 
the quality of the results. 



3. Variability 

3.1. Variability Statistics 

In order to ensure that we were using the original release of 
Kepler data correctly, we first reproduced the calculations of B 1 1 
as exactly as possible, and compared our results to theirs. We 
found no discrepancies once the different data reduction meth- 
ods had been taken into account. 

There are a variety of statistics that can be used to quantify 
variability. The Cll study use the PDC 'dispersion' which can 
be downloaded from (NStED) together with other pre-computed 
statistics, including the light curve median and reduced x^- 
Dispersion is defined as the 1 sigma rms scatter around the me- 
dian magnitude of the light curve. 
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Fig. 2. Variability index (described in Section 3J_) against Kepler magnitude for the Ra w, PD C and ARC Kepler Ql data. The solid 
line, in the same position on each graph, shows the photometric uncertainty (see Section [T2| ). The removal of true stellar variability 
at medium levels by the PDC can be seen in the dearth of stars around 7?var = 10^ - 10"^ ppm in the middle panel. 



Instead of dispersion, BIO, 11 measured the light curve 
'range', which is essentially a measure of the peak-to-peak vari- 
ation. The effect of high-frequency noise was removed either by 
smoothing the light curve on 10-hour timescales (BIO) or by dis- 
carding the upper and lower 5 percentiles (Bll). The choice of 
statistic used to study the variability is somewhat arbitrary and 
we confirm that this choice does not significantly alter the re- 
sults. 

We used the empirical three-section cut of CI 1 to distinguish 
between likely dwarfs and giants based on surface gravity log g 
and eff'ective temperature Teff , rather than the simpler log g cut 
used by BIO. It has since been noted that the KIC contains some 
misidentifications ( |Koch et al.|2010| ), but since these are not ex- 
pected to be numerous enough to aff'ect our results, and for ease 
of comparison to CI 1, we used the KIC values without modifi- 
cation. 

Based on the method of B 1 1 , we have chosen to use the range 
/?var, between the 5th and 95th percentile, for the median nor- 
malised light curve as our variability statistic. Dispersion and 
reduced provide an appropriate measure of variability that 
is believed to be primarily stochastic and Gaussian. Pulsations 
and rotational variability do not meet these criteria and there- 
fore a measurement based on the peak-to-peak variations in the 
light curve is considered more relevant. Selecting the 5th to 95th 
percentile range reduces the noise on the peak-to-peak measure- 
ment. 

/?var measurements for the Raw, PDC and ARC data are 
compared in Figures |2] and |3] The ARC clearly removes most 
systematic eff'ects, reducing the lower envelope of points to the 
photo noise limit, however it does not have the side eff'ect of sup- 
pressing intermediate amplitude variability (as done by the PDC; 
this is apparent in the scarcity of points around 7?var = 10^-10"^ 
ppm in Figure [3] and the middle panel of Figure |2]) Another un- 
fortunate side eff'ect of the PDC is the introduction of high fre- 
quency noise in some light curves, which does not occur with 
the ARC. 




10^ 
^var (ppm) 



Fig. 3. Histogram of variability for stars 13 < Kepmag < 14 for 
the Raw, PDC and ARC Kepler Ql data. 



3.2. Solar Comparison and Variability Fraction 

To compare the properties of the high and low variability stars 
we make a cut in /?var based on a comparison to twice the vari- 
ability level of the active Sun. This level is somewhat arbitrary 
but provides an appropriate way to divide stars with approxi- 
mately solar levels of variability and quieter, from those which 
are significantly more variable than the active Sun. The solar 
7?var value was calculated from the SOHO/VIRGO summed g-\-r 
light curves for the active Sun in the year 2000, because these 
provide the closest match to the Kepler bandpass (BIO). The so- 
lar 7?var value was calculated to be 766 ppm from the average 
obtained using a sliding 33 day section of data over 2 years, cen- 
tred on the activity maximum. An empirical fit to the median of 
the photometric uncertainties on the light curves in 0.5 mag bins 
provides an estimate of noise levels across the range of mag- 
nitudes. The equation of the solar equivalent line is composed 
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Fig. 4. Division of high and low variability dwarf stars (red solid 
line) at twice the solar value (red dashed line). The solar level 
is calculated from the active Sun level (cyan dashed line) and 
the noise (magenta solid line), which itself is a combination of 
background (blue dashed) and photon noise (green dashed). The 
orange dashed line marks the position of the solar line as deter- 
mined by BIO. 



of the solar 7?var measurement, photon and background noise 
terms. 



/ = solar 4ar + lO^^^^^g-^^^' + lO^-^^^^S- 



■Pif 



(3) 



where best fit values for pi and p2 are 0.4 and 8.0 respectively. 
See Figure |4] for a graphical representation of the division pro- 
cess. The orange dashed line in Figure [4]represents the solar line 
as determined by BIO, where extensive visual comparison of so- 
lar and Kepler light curves is used to estimate a solar level. 

The fractions of dwarf stars with variability greater than the 
solar level are shown in Table [T] We find that 60% of dwarf stars 
are more variable than the active Sun on 33 day timescales. BIO 
find 46% of dwarf stars to be more variable than the solar level. 
Cll define stars as significantly variable if the reduced is 
greater than 10, and find the fraction of dwarf stars meeting this 
criteria to be 18%. 

The variability fraction is strongly dependent on the choice 
and definition of the solar and photometric noise levels. A larger 
fraction of stars fall into the low variability subset when the di- 
vision of BIO is used, due to the small diff'erence in the solar 
value used and more importantly the position of the estimated 
photometric noise level. Our empirical fit to the median of the 
photometric uncertainties produces a line that more closely fol- 
lows the shape of the lower envelope in Figure |4] than the visual 
estimate of BIO. 

We calculated the random uncertainty for each of the vari- 
ability and periodicity fractions listed in Table [T] by performing 
10,000 measurements using random samples of 80% of the data, 
and found these to be < 0.5% and hence negligible in compari- 
son to diff'erence introduced by the choice of dividing line. 

It is important to note that for many stars in the low vari- 
ability sample, the lower limit of 7?var is due to the photometric 
precision of the Kepler satellite, shown as the magenta line in 
Figure [4] However, due to the position of the dividing line (solid 
red line in Figure]?]) and its dependence on the photometric pre- 
cision, there are still a large number of stars that are significantly 
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Fig. 5. Density plot of eff'ective temperature against 7?var, show- 
ing the decrease in variability with increasing temperature for all 
the selected dwarf stars (top) and for those with Kepler magni- 
tude < 14 (bottom). This illustrates that the dearth of low vari- 
ability starts at low temperatures is not a result of the increased 
noise floor of the cool, faint stars. The dashed line shows the 
solar variability level and the solid line is twice the solar level. 



more variable than the photometric noise limit for their magni- 
tude. There are -1,500 stars in the low variability sample that 
are above twice the photometric noise level. 



3.3. Comparison of Low and High Variability Stellar 
Properties 

BIO, 11 and Cll highlighted various relationships between vari- 
ability and stellar properties. We revisit these using the ARC 
data. As CI 1 noted, there is a decrease in variability with temper- 
ature which cannot be entirely explained by the increased noise 
levels in the fainter cool stars (see Figure [5]). Pearson's correla- 
tion coefficient of eff'ective temperature and logiQ(7?var) for the 
whole dwarf sample is -0.31. 

To examine whether the high and low variability samples be- 
long to diff'erent stellar populations, we investigated the possi- 
ble diff'erences in their spatial distribution and kinematics. Due 
to dynamical heating of the galactic disk over long timescales, 
older populations of stars tend to have larger galactic scale 
heights and l arger velocity dispersions than yo unger popula- 
tions (see e.g. [Freeman & Bland-Hawthorn|2002[ and references 
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Table 1. Fractions of stars m ore v ariable than once and twice the solar level, and fractions of periodic stars in the high and low 
variability groups (see Section 4.1 for definition). The random uncertainty for each of the measurements is < 0.5%, and is negligible 
in comparison to diff'erence introduced by the choice of dividing line, demonstrated by the 1 and 2 times the solar value comparisons 
listed in this table (see Section[T2|for further discussion). 
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Fig. 6. Typical light curves from the high and low variability groups {left and right respectively) for each spectral type. 



therein). We therefore examined the distribution of the two sam- 
ples in galactic coordinates, which is shown in Figure |7] 

The division in /?var introduces a magnitude bias between 
the high and low samples, with a greater proportion of the low 
variability objects at higher magnitudes. To reduce any possible 
eff'ect of this bias on our comparison tests, we selected high and 
low variability subsets with approximately equal magnitude dis- 
tributions. This was done by dividing the stars into 1 mag bins 
(with the exception of stars brighter than 8th and fainter than 
16th in Kepler magnitude, which were treated separately), and 
choosing a random set from the more populated bin to match 
the number in the smaller one. This process, used for the spatial 
distribution and proper motion comparison tests, also serves to 
select an equal number of stars in the high and low variability 
sets. 

When using the PDC data, a diff'erence in spatial distribution 
of the two subsets is evident, as noted by CI 1. The low variabil- 
ity stars appear approximately uniformly distributed, whereas 
the high variability stars are concentrated towards low galactic 



latitudes. Using the PDC data alone, one could hypothesise that 
the variable sample may correspond to a younger, thin disk pop- 
ulation, and the quiet sample to an older population with larger 
scale height. Cll hinted at this conclusion but argued that this 
could also be an artefact of increased crowding at low galactic 
latitudes, causing higher levels of photometric dispersion. 

By visually comparing the spatial distribution of the high 
and low variability samples from the PDC and ARC, it becomes 
evident that the eff'ect is dependant on the choice of reduction 
method (Figure[7|. This suggests that contamination eff'ects were 
the dominant cause of the apparent higher variability close to the 
galactic plane, and that the ARC is better at removing the asso- 
ciated systematics. There is still some sign that active F stars 
are more concentrated at lower latitudes (Figure [8]), but it is too 
subtle to draw any strong conclusions. 

The giant contamination in the F stars is higher than other 
spectral classes because they span the temperature range where 
the giant branch intersects the main sequence on the HR dia- 
gram. This may lead to the apparent concentration of variable F 
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Fig. 7. The spatial distribution of low (left) and high (right) vari- 
ability stars for ARC (top) and PDC (bottom). Each has been 
subject to a random selection of equal numbers of stars per mag- 
nitude bin, before a random selection of 6000 for each panel was 
selected. The PDC panels show a slight tendency for the high 
variability sample to be more concentrated towards the galactic 
plane, whereas this effect is not significant in the ARC, suggest- 
ing it is more eff'ective at removing contamination related sys- 
tematics. 



stars at low galactic latitudes, if the variable and luminous gi- 
ants are visible out to greater distances in the galactic plane, and 
hence appear with higher density in these regions. 

We also performed a more quantitative two- sample 
Kolmogorov-Smirnov (KS) test. The results are shown in 
Table |5| and confirm that the galactic latitudes of the low and 
high variability F stars are very unlikely to be drawn from the 
same distribution, compared to a high probability for the M 
stars. The distribution shape parameters, also shown in Table [2j 
display the trends of increasing median galactic latitude with 
later spectral type and the subtle diff'erences between high and 
low variability subsets, although these are not sufficient to prove 
the stars belong to diff'erent populations. The spatial distribution 
plots also serve to show that the correction applied using ARC 
aff'ects all areas of the CCD equally, despite being implemented 
on each mod.out separately. 

To examine the eff'ect of crowding, we used the contamina- 
tion fraction available from the Kepler mission archive, which 
provides an approximate measure of the fraction of flux at- 
tributed to the target in a 21 x 21 pixel aperture. Values range 
from 0, implying no contamination, to 1, which indicates the 
flux is essentially all background. There is a weak correlation be- 
tween contamination fraction and 7?var, shown in Figure [9] We 
tested whether placing a constraint on the contamination frac- 
tion, using only targets in the range 0.1 to 0.2, altered the ap- 
pearance of Figure [T] but it did not. 

A potential indication that the low and high variability sub- 
sets belong to difl'erent stellar populations can be seen in the 
proper motion distributions (Figure [TO]), for stars with total 
proper motion greater than zero. It shows a weak but noteworthy 
trend for the high variability group to have lower proper motion 
than the low variability group. We also tested the significance of 
these diff'erences using two-sample KS tests, and parameterised 
the shape of the distributions, the results of which can be seen in 
Table [3] The KS test results show that with the exception of A 




10 11 12 13 14 15 16 17 18 

h (deg) 

Fig. 8. Histogram showing the galactic latitude distribution of 
low (solid line) and high (dashed line) variability stars for each 
spectral type. Sample selection ensured the orientation of the 
Kepler field on the plane did not introduce biases. An equal sam- 
ple of high and low variability stars from each magnitude bin 
was also selected. The increased number of high variability F 
stars at low galactic latitudes may arise from giant contamina- 
tion of the sample (see Section [33]). 



stars, and to some extent F, the proper motion values for the high 
and low variability groups are very unlikely to be drawn from 
the same distribution, while the shape parameters highlight more 
subtle diff'erences in shape and trends between spectral types. 
This is consistent with the view that higher variability stars are 
younger and therefore have lower proper motions. This conclu- 
sion is dependant on the stars being at the same distance and lo- 
cation on the sky. In our samples, the magnitude and spatial dis- 
tributions are approximately equal for the high and low variabil- 
ity subsets, which should satisfy this condition. One caveat ap- 
plicable to the proper motion distributions is that the giant con- 
tamination in the M dwarfs (described by CIO) could increase 
the apparent number of high variability stars with low proper 
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Table 2. Statistics of galactic latitude (b) distributions, including two-sample Kolmogorov-Smirnov tests, median, median absolute 
deviation (MAD) and skew of the low (L) and high (H) variability samples. 



Sample KS MedL/? MedH/? MAD L/? MAD m Skew Skew HZ? 



A 


0.016 


11.68 


12.00 


1.76 


1.57 


0.70 


0.54 


F 


1.4x10-9 


12.52 


11.79 


1.92 


1.70 


0.35 


0.59 


G 


0.002 


12.88 


12.65 


2.08 


1.95 


0.21 


0.28 


K 


0.1 


13.29 


13.32 


2.21 


2.08 


0.07 


0.04 


M 


1.0 


13.65 


13.58 


2.21 


2.08 


-0.004 


-0.016 



Table 3. Statistics of proper motion distributions, including two-sample Kolmogorov-Smirnov tests, median, median absolute devi- 
ation (MAD) and skew of the low (L) and high (H) variability samples. 



Sample /^^total M^d l^pm^^^^i Med Yipm^^^^i MAD l^pm^^^^i MAD Yipm^^^^i Skew l^pm^^^^i Skew H/7m^Q^^| 



A 0.61 0.0057 0.0063 0.0032 0.0030 4.61 1.57 

F 0.09 0.0077 0.0076 0.0038 0.0033 9.57 12.66 

G 0.0009 0.010 0.0091 0.0045 0.0037 12.06 4.45 

K 0.002 0.013 0.011 0.0061 0.0045 8.75 9.41 

M 1.5x10-6 0.033 0.020 0.013 0.0091 2.41 4.48 




10- 



10^ 10^ 
^var (ppm) 



10^ 



Fig. 9. A weak correlation can be seen between contamination 
fraction and /?var, for non-zero contamination values. 



motion. We have also considered potential effects of Malmquist 
bias in our results, but believe that by selecting equal magnitude 
distributions and comparing each spectral type separately, any 
effect introduced should be negligible. 



4. Periodic and Stochastic Variability 

We investigated the nature of the variability, comparing the low 
and high variability samples, using a combination of visual ex- 
amination, the statistics described in Section|3j and periodogram 
analysis. The typical characteristics of the light curves also vary 
between spectral types as shown in the examples in Figure [5] 
Visual examination of the light curves shows that A and F stars 
contain many pulsators with a large range of amplitudes, while 
later types show rotational modulation and more stochastic vari- 
ability. This trend is also apparent in the periodicity and stochas- 
ticity tests described in this section. 



4.1. Periodicity 

For each light curve, we computed a periodogram by least- 
squares fitting of a sinusoid plus a constant at each trial period. 
The periodogram is expressed in terms of the statistic 



{xl- 



■X' 



(4) 



where xl is the reduced chi-squared of the light curve with re- 
spect to a constant value, and;^^ is the reduced chi-squared with 
respect to the best-fit sinusoid. We used 500 logarithmically- 
spaced periods between 0.01 and 100 days. However, only pe- 
riods between 1 hour (twice the sampling rate) and 16 days (half 
the time span of the data) were considered valid. We chose this 
approach over limiting the initial search space to this range, be- 
cause this would have led to many spurious detections close to 
the upper limit. Note that some objects with true periods above 
16 days may be detected at harmonics of their true period. 

The dearth of stars with periods detected close to the 16 day 
limit is a bias introduced by the stringent period identification 
method. Many objects that would occupy this period range have 
broad periodogram peaks at around the 16 day limit and the se- 
lection method has been designed to only accept periodogram 
peaks that drop to the median level before reaching the edges of 
the allowed range. This prevents false detections of periods at 
the length of the permitted dataset, and misidentification of the 
highe st periodogram continuum point as a genuine peak. 

Bretthorst'('1998) describes the Bayesian theory of modelling 
data using a sinusoid plus white noise. In this formalism, our 
metric can be shown to be 



S =2C(oj)/{Ncrj^^), 
where(^^) = [R(ajf -h I(ojf) /N, 



(5) 



(6) 



is the squared magnitude of the discrete Fourier transform of the 
data. N is the number of data points in the light curve, R(oj) = 
Yjf^i cos(ojti), = 2/1 1 sin(6L)^/) and crrnis is the rms scatter 
of the data. We have assumed the data is mean- subtracted. 

As discussed in B 1 1 , the choice of a threshold for periodic- 
ity detection is best made empirically. We selected S = 0.3 as 
an appropriate threshold between clear-cut periodic variability 
and non-periodic or ambiguous cases. This threshold is slightly 
lower than the value of S = 0.4 used in the Monitor project (see 
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Proper Motion Histograms 




0.020 



0.04 



0.05 



0.05 



0.02 0.04 0.06 0,08 0.10 0.12 
Total Proper Motion (7yr) 

Fig. 10. Distribution of proper motion values for the low (solid 
line) and high (dashed line) variability samples, where proper 
motion > 0, in a box of even galactic latitude and longitude dis- 
tribution, with equal numbers of high and low variability stars 
selected in each magnitude bin. The uncertainty on the proper 
motion values is 0.02"/year but given the high numbers in each 
sample these results are still significant. 



e.g. [Irwin et al.|2QQ6| , which is reasonable given the vastly su- 
perior time-sampling of the Kepler data. 

This threshold can be compared to a power spectrum cut, 
since the power spectrum is given by P S(6l)) = 4C(co)/N where, 
following [Kjeldsen & Bedding] ( |1995| ) we have normalised the 
power spectra such that a sinusoidal oscillation of amplitude A 
gives rise to a peak of height in the power spectrum. Our 
threshold is then equivalent to a power spectrum threshold of 



^PS ~ ^-^^rms- comparison Kjeldsen & Bedding ( 1995| 



01 



'rms- ■ 

give an expression for the noise level in the power spectrum 

(Tp3 = 4(Trms/A^- 

In our data, A/^ » 1, so our threshold is conservative when 
comparing to white noise. However, white noise is not the dom- 
inant factor defining our ability to detect periodicities. To test 



the appropriateness of our periodicity threshold on data with re- 
alistic noise properties, we ran a set of simulations where we 
injected periodic signals into actual Kepler light curves. We 
randomly selected 1000 Ql light curves with low variability 
(7?var < ^var,Sun) and no significant period (S < 0.25), and in- 
jected sinusoidal signals into them, with random periods uni- 
formly distributed between 2 and 16 days, and random ampli- 
tudes ranging from 0.1 to 10 times the high-frequency noise level 
(measured as the scatter in the diff'erence between consecutive 
flux measurements), with a distribution of amplitudes relative to 
the noise level that was uniform in log. We then stored the best- 
fit value of S and the corresponding period. We found that S was 
always < 0.3 when the best-fit period diff'ered by > 10% from the 
injected value, suggesting a false detection rate < 10"^. As one 
would expect, the missed detection rate is period- and amplitude- 
dependent: there were only 5 cases in our simulations where the 
recovered period was within 10% of the injected one but S was 
< 0.3, and these were all for /^injected > 12 days and amplitudes 
significantly smaller than the noise. However, our simulations 
did not include a large enough number of these low- amplitude, 
difficult to detect signals to make a more detailed estimate of the 
completeness. Since the period search in Ql data was intended 
as a preliminary exercise, and the results will become much more 
reliable when additional data is included, we decided additional 
simulations were beyond the scope of the present paper. 

Using the criteria outlined here, 16% of the dwarf stars 
are determined to be periodic, and the fraction for each spec- 
tral type is given in Table [T] We note that these numbers vary 
slightly from those quoted in Bll and Cll. This indicates that 
our threshold for periodicity is more stringent, designed to de- 
tect only genuine periodicities (or strong harmonics) within the 
limits of the frequency resolution. The periodic fraction will un- 
doubtably increase with longer datasets and more efficient quasi- 
periodic signal detection methods. 

Histograms of the detected periods for each spectral class 
are shown Figure [TT] When interpreting these, one should bear 
in mind that our period sensitivity is non-uniform, and that the 
largest peak in the periodogram is not necessarily at the true 
period but can be at one of its harmonics. The distribution is 
also truncated close to the maximum period due to the defini- 
tion of 'peak' required for selection. These period distributions 
should thus be taken as indicative only (the period sensitivity 
will improve vastly with longer time coverage). Nonetheless, 
differences between the histograms are obvious, with the typi- 
cal period clearly increasing towards later spectral types. 

The majority of A stars, and about half of F stars, have very 
short periods (< 2 days). These are likely to be pulsators (many 
of the A and F stars are located in the instability strip), although 
they could also be unpublished close binaries (where ellipsoidal 
variations and mutual heating induce sinusoidal variations) or 
active stars with very short rotation periods. On the other hand, 
G, K and M stars, along with the rest of the F stars, have sig- 
nificantly longer periods (> 5 days), as one might expect from 
rotational modulation of active regions. This distinction is con- 
firmed by the appearance of the light curves (Figure |6]). The 
trend of longer periods towards later spectral types is also vis- 
ible amongst the late type objects. 

4.2. Degree of Stochasticity 

A significant difference in the nature of variability across diff'er- 
ent spectral types is shown by the degree of stochasticity. We per- 
formed a simple measure of this based on the number of peaks, 
A^pk, in the sine-fitting periodogram which have power greater 
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Fig. 11. Period distribution for each spectral type, for stars where 
a significant period between Ihr and 16d has been detected (see 
Section [4rT] for method and caveats), showing an increase in pe- 
riod towards later spectral types. 



than 10% the maximum power. We selected only stars above the 
periodicity threshold described in Section 4.1 An example A/pk 
determination is shown in Figure 12 



The histograms showing the distribution for each spec- 
tral type are displayed in Figure [13] The light curves of pulsat- 
ing stars are dominated by near-sinusoidal variability at one, or 
a few, clear dominant periodicities, corresponding to low val- 
ues of A^pk. Rotational variables are less well modelled by sine- 
fitting. They typically have a more complex periodogram, with 
significant power at multiple harmonics of the dominant pe- 
riod (see e.g. |Boisse et al.|2 009 ) and some stochastic variability 
(with p ower at all frequenci es, as seen in the case of the Sun, 
see e.g. |Aigrain et al.||20Q4| ). This results in somewhat larger 
values of A^pk. The high variability group show a tendency to- 
wards slightly lower levels of stochasticity than the low variabil- 
ity group, which may arise from the strongly periodic pulsating 
stars that typically have high amplitude variations. It is impor- 
tant to note that the A^pk statistic is not a quantitive measure of 
stochasticity, but intended for use as a comparison between the 
low and high dispersion samples. 

Alternative metrics for measuring timescale, periodicity and 



stochasticity were demonstrated by |Walkowicz & Basri| ( |201Q 



who use time separation between points where the diff'eren- 
tial light curve crosses zero, and the time separation between 
changes in the sign of the slope in the light curve. This method, 
without smoothing of the light curve, is slightly more sensitive 
to noise than the A/^pk approach. 

4.3. Harvey Model Fitting 

We also studied the typical stochastic properties of each spec- 
tral class by fitting a model to the average power spectra of 
each. The Fast Fourier Transform (FFT) of each light curve was 
computed and the median, 10th and 90th percentile at each fre- 
quency were used to create the plots in Figure 14 These spectra 
were smoothed using a nonlinear median boxcar filter to remove 
the effects introduced by the variable pointing stars (evidence of 
which can be seen in the peaks of the cyan line in Figure [14]). 
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Fig. 13. Distribution of the number of peaks with power greater 
than 10% of the maximum peak for low (solid line) and high 
dashed line) variability samples, where S > 0.3. 

Autoregressive (AR) models are commonly used to describe 
stochastic processes and take the form. 



(7) 



where Xt and Xt-i are data values at time index t and t-1 re- 
spectively, c is a constant, p is the order of the model, cp are the 
parameters and e is white noise. Such a process has a spectral 
density 



(Btvy 



(8) 



where v is frequency. A/ is the amplitude of the /th compo- 
nent, Bi is its characteristic timescale, and C/ is the slope of the 
power law. By considering the power spectra as a sum of A^ dif- 
ferent AR(1) spectra, one can fit a model of this form, first in- 
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Fig. 12. Graphical representation of the calculation of A^pk for periodic {left) and non-periodic (right) stars. The top panels show the 
original power spectrum, smoothed spectrum, and the maximum peak (dashed line). The original power spectrum is divided by the 
smoothed one to get the corrected spectrum in the centre panels. Peaks above 10% of the maximum value (horizontal dashed line) 
are then counted. The bottom panels shows the light curve (solid line) with best fit sinusoid (dashed). 



troduced by 'Harvey" ("1985) and previously used to model stellar 
spectra (see e.g. Aigrain et al.,2QQ4{pichel et al.|2Q09l ). 

The Harvey model consists of a constant background level, a 
component to describe the intrinsic power at each frequency, and 
a second component to describe the extra power at ~ 8 x 10"^ Hz 
introduced by the onboard motor vibrations. The timescale and 
slope of the motor component have been set to the mean value 
found when the parameters are allowed to vary, because there 
is no reason they should be diff'erent between spectral types. 
The amplitude for each is allowed to vary because it will have a 
greater eff'ect for the fainter stars. 

The median power spectrum and fitted powerlaw models de- 
scribe the typical stochastic background variability of the stars. 
Clear trends through the spectral classes can be seen in the fit 
parameters, listed in Table [4] As expected, the background con- 
stant increases steadily towards the fainter classes. The ampli- 
tude, timescale and slope of the powerlaw fit also increases to- 
wards later spectral types. This may arise from the slower typical 
rotation periods associated with later type stars, creating large, 
slowly evolving active regions. The timescale of the A stars does 
not fit this trend, which is most likely due to the large number of 
pulsating stars shifting the median power towards slightly longer 
timescales, and the fact that the frequency resolution is not suf- 
ficient to resolve the powerlaw turnoff'. 



5. Discussion and Conclusions 

We have developed an astrophysically robust systematics correc- 
tion method (ARC), which has been shown to remove the signif- 
icant systematic trends from the Kepler Ql light curves, while 
maintaining intrinsic stellar variability. Figure [2] shows that the 
ARC performs to the same standard as the PDC in reducing the 
lower envelope of the measured variability down to the photo- 



Table 4. Best fit parameters for the Harvey models shown in 



Figure 14 
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Al 


Bi 
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A2 
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C2 


A 


L9e7 
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1.0 
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1.6e3 


4.8 


F 
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L2e5 


1.7 


9.1e6 
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4.8 
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L2e8 


l.lelO 
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2.1 
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1.6e3 


4.8 
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2.0e8 


1.2ell 
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4.2e7 


1.6e3 


4.8 


M 


3.2e8 


L9el2 


3.6e5 


2.8 


1.2e8 


1.6e3 


4.8 



metric uncertainty. Crucially, unlike the PDC it does not remove 
real astrophysical signatures (predominantly at medium variabil- 
ity levels), leading to the apparent bimodality in variability level 
that can be seen in the PDC panel of Figure [2] The correction is 
applied individually to each mod.out but the results are shown to 
be consistent across the whole field of view (Figure [7]). 

Following Bll, we quantified the variability level of each 
star as the interval between the 5* and 95* percentile of its light 
curve, and we focused our study on the main-sequence stars, 
separating them from giants using the Teff -dependent log g cut 
defined by CIO. A division between low and high variability 
was made using a magnitude dependent cut in ^var at twice 
the equivalent variability measurement for the active Sun, re- 
sulting in a split of 80% low and 20% high. We examined the 
relationship between the nature of the variability and physical 
properties of the stars, and compared our work to that of Cll, 
who use the PDC data and B 10,1 1, who use a polynomial fitting 
correction method. As suggested by BIO and Cll, we see clear 
evidence for a decrease in variability with increasing tempera- 
ture (Figure |5]), which cannot be entirely explained by selection 
eff'ects or the noise floor at fainter magnitudes. This is equiva- 
lent to an increase of variability with increasing 5 - V, as noted 
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Fig. 14. Harvey model fits to the median periodogram for each spectral type. The blue shaded area represents the region between 
10-90th percentiles of variability. The cyan line shows the unsmoothed median spectrum, and the fit (red dashed line) is made to the 
smoothed spectra (black line). The individual components of the fit are shown as blue dashed lines (see text for full explanation). 



by Isaacson & Fischer (2010), who use S^^^i indicator of 
chromosphenc activity. 

The spatial distribution and proper motion of the stars were 
used to look for an indication that the low and high variabil- 
ity samples come from diff'erent stellar populations. Cll sug- 
gest that the high variability sample belongs to a younger pop- 
ulation, existing closer to the galactic plane than the older low 
variability sample. We note a weak correlation between contam- 
ination fraction and variability, which could potentially increase 
the apparent variability in the more contaminated regions close 
to the galactic plane. The proper motion distributions displayed 
in Figure 10 indicate that, for the M stars at least, more ac- 
tive stars tend to have slightly lower proper motions. This may 
suggest they are younger, although distance estimates and radial 
velocity measurements would be needed to interpret the proper 
motion distributions robustly, and it should be noted that giant 
contamination could aff'ect this result. 

There may exist interesting further correlations between the 
variability level and other parameters listed in the KIC. For ex- 
ample, we compared the surface gravity and metallicity at dif- 
ferent variability levels, and did find apparently statistical sig- 
nificant diff'erences (in K-S test terms). However, we chose not 
to report them here, because these parameters have larger uncer- 
tainties than Teff , and are eff'ected by non-negligible biases as a 
function of Galactic positi on (because of crowding and extinc- 
tion, see e.g. Vemer et"aL][2011 , and references therein), which 
cannot easily be disentangled from actual diff'erences between 
populations. Any biases aflfecting the \ogg information in the 



KIC may also have an indirect eff'ect on the present study be- 
cause they_£ouMaflfec^he^^ dwarfs and gi- 
ants. |^mereFar](|20TT); Brown et al. (2011 ) investigated the 
reliability of the KIC parameters in detail, and concluded that 
the Teff and log g estimates are reliable within the stated uncer- 
tainties. This implies that any biases should have only a very 
minor eff'ect on the Cll dwarf/giant cut which we used. 

One of the most interesting results in this work is the pe- 
riod distribution for each spectral type, displayed in Figure [TT] 
Although this should be considered a preliminary test due to 
the short timespan of the dataset, there is clear evidence that, 
for light curves with a clear periodic component, the period in- 
creases towards later spectral types. The fraction of periodic 
stars also varies with spectral class, increasing towards later 
types to a maximum in the K stars, before reducing again in 
M stars. The relatively high number of stars in which we were 
able to identify a periodic or quasi periodic modulation in this 
dataset alone suggests that, when further quarters of Kepler data 
are added, it should be possible to measure periods for a signifi- 
cant fraction of all Kepler stars - and thus calibrate the evolution 
of angular momentum for intermediate and low-mass stars on the 
main sequence to an unprecedented level. 

We also investigated the stochastic component of the vari- 
ability in two diff'erent ways. First, we measured the number of 
periodogram peaks with power greater than 10% of the maxi- 
mum for the stars which pass our periodicity selection threshold. 
For all spectral type, the distribution of this statistic, which we 
call A^pk, has a peak at low values, corresponding to clearly pe- 
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riodic light curves dominated by a single frequency, or a small 
number of frequencies, as would be expected for pulsating stars. 
However, there is also another peak at high A^pk, whose ampli- 
tude increases towards later spectral types, reaching a maximum 
in the K stars and decreasing slightly in the M stars (as did the 
periodicity fraction). This peak corresponds to quasi-periodic 
light curves, as expected for rotational variables with evolving 
active regions. 

We also parameterised the stochastic component of variabil- 
ity was performed by fitting auto-regressive models (also known 
as Harvey models) to the median power spectrum for each spec- 
tral type. Because of the relatively limited frequency resolution 
(caused by the short duration of the dataset used in this study), 
the fitted model parameters cannot be considered definitive, but 
they enable a preliminary inter-comparison. We find that the typ- 
ical amplitude, timescale and power-law index all increase to- 
wards later types. 

This is consistent with our other tests, and supports a broadly 
coherent picture of main- sequence variability. The hotter, ear- 
lier spectral classes show a lower variability level on the whole, 
and the variables tend to show clearly periodic behaviour on 
short time- scales, as expected from pulsations. The shape of 
their power spectra suggest that these stars possess smaller ac- 
tive regions that evolve more quickly. By contrast, the cooler, 
later type stars show larger amplitude variability on longer time- 
scales, with quasi-periodic rather than periodic behaviour, and 
appear to possess more slowly evolving, larger active regions. 
The K stars have the highest fraction of light curves with a sig- 
nificant periodicity, but these are also the most complex, with the 
largest number of significant frequencies per object. 

We are currently working towards a systematics correction 
suitable for the Q2 data, which requires treatment for discon- 
tinuities and safe mode eff'ects in conjunction with the ARC 
method. Once complete, this will allow us to perform a simi- 
lar investigation on a much longer dataset, providing vastly im- 
proved constraints on the periodicity measurements, and reveal- 
ing changes in the variability though time (i.e. pulsation and ac- 
tivity cycles). These measurements will improve with each new 
release of Kepler data. We will also seek to compare our vari- 
ability measurements to other metrics and timescales, for exam- 



ple the CDPP on transit timescales calculated by [Gilliland et al. 
( 2Q11| ). Finally, the results presented in this paper also have im- 
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